library(survey)
library(dplyr)
library(ggplot2)
library(gridExtra)

# Chapter 10

### Figure 10.1
summary(factor(df$Q21)) # Iran
summary(factor(df$Q25)) # Transgender ban
summary(factor(df$Q11)) 

df2 <- df
df2 <- df2[df2$Q11 < 77,]

# Attack Iran treatments
df2$P_ASSIGN2A <- NA
df2$P_ASSIGN2A[df2$P_ASSIGN2==1] <- 1 # Control
df2$P_ASSIGN2A[df2$P_ASSIGN2==5] <- 2 # Oppose
df2$P_ASSIGN2A[df2$P_ASSIGN2==4] <- 3 # Support

# Transgender Ban treatments
# P_ASSIGN3 == 1    # Control
# P_ASSIGN3 == 2    # Oppose
# P_ASSIGN3 == 3    # Support
# P_ASSIGN3 == 4    # Divided

# Q21T
df2$Q21T <- NA
df2$Q21T[df2$Q21 < 3] <- 1
df2$Q21T[df2$Q21 > 2 & df2$Q21 < 77] <- 0

# Q25T
df2$Q25T <- NA
df2$Q25T[df2$Q25 < 3] <- 1
df2$Q25T[df2$Q25 > 2 & df2$Q25 < 77] <- 0

# Create binaries for conditions
df2$isupp <- NA
df2$isupp[df2$P_ASSIGN2A==1] <- 0
df2$isupp[df2$P_ASSIGN2A==3] <- 1

df2$iopp <- NA
df2$iopp[df2$P_ASSIGN2A==1] <- 0
df2$iopp[df2$P_ASSIGN2A==2] <- 1

df2$tbo <- NA
df2$tbo[df2$P_ASSIGN3==1] <- 0
df2$tbo[df2$P_ASSIGN3==2] <- 1

df2$tbs <- NA
df2$tbs[df2$P_ASSIGN3==1] <- 0
df2$tbs[df2$P_ASSIGN3==3] <- 1

df2$tbd <- NA
df2$tbd[df2$P_ASSIGN3==1] <- 0
df2$tbd[df2$P_ASSIGN3==4] <- 1


# Survey designs
w2_design <-
  svydesign(
    id = ~ 1,
    weights = ~ weight2,
    data = df2
  )

# T-tests
# Iran
iopp <- svyttest(Q21T~iopp, w2_design, na.rm = TRUE)
isup <- svyttest(Q21T~isupp, w2_design, na.rm = TRUE)

# Transgender ban
tbo <- svyttest(Q25T~tbo, w2_design, na.rm = TRUE)
tbs <- svyttest(Q25T~tbs, w2_design, na.rm = TRUE)
tbd <- svyttest(Q25T~tbd, w2_design, na.rm = TRUE)

# Build out plots
# Iran
pt.est <- c(isup$estimate, 
            iopp$estimate)
ci.low <- c(isup$conf.int[1], 
            iopp$conf.int[1])
ci.high <- c(isup$conf.int[2], 
             iopp$conf.int[2])
plot.labs <- c("Support","Oppose")
df.plot <- data.frame(Treatment=plot.labs, estimate=pt.est, 
                      lower=ci.low, upper=ci.high)
df.plot$Treatment <- factor(df.plot$Treatment, levels = plot.labs)

g.i <- ggplot()+
  geom_pointrange(data=df.plot, 
                  mapping = aes(x=Treatment,
                                y=estimate,
                                ymin=lower,
                                ymax=upper),
                  size = 0.5, color = "black")+
  geom_hline(yintercept = 0, color = "red")+
  ylim(-.2,.2)+
  ylab("Difference in Means from Control Condition")+
  xlab("Treatment Condition")+
  ggtitle("Iran Strikes")+
  theme_bw()

# TB
pt.est <- c(tbs$estimate, 
            tbo$estimate, 
            tbd$estimate)
ci.low <- c(tbs$conf.int[1], 
            tbo$conf.int[1],
            tbd$conf.int[1])
ci.high <- c(tbs$conf.int[2], 
             tbo$conf.int[2],
             tbd$conf.int[2])
plot.labs <- c("Support","Oppose", "Divided")
df.plot <- data.frame(Treatment=plot.labs, estimate=pt.est, 
                      lower=ci.low, upper=ci.high)
df.plot$Treatment <- factor(df.plot$Treatment, levels = plot.labs)

g.tb <- ggplot()+
  geom_pointrange(data=df.plot, 
                  mapping = aes(x=Treatment,
                                y=estimate,
                                ymin=lower,
                                ymax=upper),
                  size = 0.5, color = "black")+
  geom_hline(yintercept = 0, color = "red")+
  ylim(-.2,.2)+
  ylab("Difference in Means from Control Condition")+
  xlab("Treatment Condition")+
  ggtitle("Transgender Ban")+
  theme_bw()

grid.arrange(g.i,g.tb,ncol=2)



############################################################################



### Figure 10.2
summary(factor(df$Q27))
summary(factor(df$Q11)) 

df2 <- df
df2 <- df2[df2$Q11 < 77,]

# Q27T
df2$Q27T <- NA
df2$Q27T[df2$Q27 < 3] <- 1
df2$Q27T[df2$Q27 > 2 & df2$Q27 < 77] <- 0

# Guns in schools P_ASSIGN4
# 1 - Control
# 2 - More Safe
# 3 - Less Safe
# 4 - Divided

df2$gms <- NA
df2$gms[df2$P_ASSIGN4==1] <- 0
df2$gms[df2$P_ASSIGN4==2] <- 1

df2$gls <- NA
df2$gls[df2$P_ASSIGN4==1] <- 0
df2$gls[df2$P_ASSIGN4==3] <- 1

df2$gd <- NA
df2$gd[df2$P_ASSIGN4==1] <- 0
df2$gd[df2$P_ASSIGN4==4] <- 1


# Survey designs
w2_design <-
  svydesign(
    id = ~ 1,
    weights = ~ weight2,
    data = df2
  )


# T-tests
# guns in schools
gms <- svyttest(Q27T~gms, w2_design, na.rm = TRUE)
gls <- svyttest(Q27T~gls, w2_design, na.rm = TRUE)
gd <- svyttest(Q27T~gd, w2_design, na.rm = TRUE)

# Build out plots
# Iran
pt.est <- c(gms$estimate, 
            gls$estimate,
            gd$estimate)
ci.low <- c(gms$conf.int[1], 
            gls$conf.int[1],
            gd$conf.int[1])
ci.high <- c(gms$conf.int[2], 
             gls$conf.int[2],
             gd$conf.int[2])
plot.labs <- c("More Safe","Less Safe", "Divided")
df.plot <- data.frame(Treatment=plot.labs, estimate=pt.est, 
                      lower=ci.low, upper=ci.high)
df.plot$Treatment <- factor(df.plot$Treatment, levels = plot.labs)

ggplot()+
  geom_pointrange(data=df.plot, 
                  mapping = aes(x=Treatment,
                                y=estimate,
                                ymin=lower,
                                ymax=upper),
                  size = 0.5, color = "black")+
  geom_hline(yintercept = 0, color = "red")+
  ylim(-.1,.1)+
  ylab("Difference in Means from Control Condition")+
  xlab("Treatment Condition")+
  theme_bw()



############################################################################



### Figure 10.3
summary(factor(df$Q21)) # Iran
summary(factor(df$Q25)) # Transgender ban
summary(factor(df$Q11T)) 

df2 <- df
df2 <- df2[df2$Q11 < 77,]

# Attack Iran treatments
# P_ASSIGN2 == 1 - Control
# P_ASSIGN2 == 4 - Support
# P_ASSIGN2 == 5 - Oppose

# Transgender Ban treatments
# P_ASSIGN3 == 1    # Control
# P_ASSIGN3 == 3    # Support
# P_ASSIGN3 == 2    # Oppose
# P_ASSIGN3 == 4    # Divided

# Q21T
df2$Q21T <- NA
df2$Q21T[df2$Q21 < 3] <- 1
df2$Q21T[df2$Q21 > 2 & df2$Q21 < 77] <- 0

# Q25T
df2$Q25T <- NA
df2$Q25T[df2$Q25 < 3] <- 1
df2$Q25T[df2$Q25 > 2 & df2$Q25 < 77] <- 0

# Survey designs
w2_design <-
  svydesign(
    id = ~ 1,
    weights = ~ weight2,
    data = df2
  )


# Iran
iran <- svyby(~Q21T,~P_ASSIGN2+Q11T, w2_design, svymean, na.rm = TRUE)

# Transgender ban
tban <- svyby(~Q25T,~P_ASSIGN3+Q11T, w2_design, svymean, na.rm = TRUE)

# Build out plots
# Iran
pt.est <- c(iran$Q21T[iran$P_ASSIGN2==1& iran$Q11T==0], 
            iran$Q21T[iran$P_ASSIGN2==4& iran$Q11T==0],
            iran$Q21T[iran$P_ASSIGN2==5& iran$Q11T==0],
            0,
            iran$Q21T[iran$P_ASSIGN2==1& iran$Q11T==1],
            iran$Q21T[iran$P_ASSIGN2==4& iran$Q11T==1],
            iran$Q21T[iran$P_ASSIGN2==5& iran$Q11T==1],
            tban$Q25T[tban$P_ASSIGN3==1& tban$Q11T==0], 
            tban$Q25T[tban$P_ASSIGN3==3& tban$Q11T==0],
            tban$Q25T[tban$P_ASSIGN3==2& tban$Q11T==0],
            tban$Q25T[tban$P_ASSIGN3==4& tban$Q11T==0],
            0,
            tban$Q25T[tban$P_ASSIGN3==1& tban$Q11T==1], 
            tban$Q25T[tban$P_ASSIGN3==3& tban$Q11T==1],
            tban$Q25T[tban$P_ASSIGN3==2& tban$Q11T==1],
            tban$Q25T[tban$P_ASSIGN3==4& tban$Q11T==1])
ci.low <- c(iran$Q21T[iran$P_ASSIGN2==1& iran$Q11T==0]-1.96*iran$se[iran$P_ASSIGN2==1& iran$Q11T==0], 
            iran$Q21T[iran$P_ASSIGN2==4& iran$Q11T==0]-1.96*iran$se[iran$P_ASSIGN2==4& iran$Q11T==0],
            iran$Q21T[iran$P_ASSIGN2==5& iran$Q11T==0]-1.96*iran$se[iran$P_ASSIGN2==5& iran$Q11T==0],0,
            iran$Q21T[iran$P_ASSIGN2==1& iran$Q11T==1]-1.96*iran$se[iran$P_ASSIGN2==1& iran$Q11T==1],
            iran$Q21T[iran$P_ASSIGN2==4& iran$Q11T==1]-1.96*iran$se[iran$P_ASSIGN2==4& iran$Q11T==1],
            iran$Q21T[iran$P_ASSIGN2==5& iran$Q11T==1]-1.96*iran$se[iran$P_ASSIGN2==5& iran$Q11T==1],
            tban$Q25T[tban$P_ASSIGN3==1& tban$Q11T==0]-1.96*tban$se[tban$P_ASSIGN3==1& tban$Q11T==0], 
            tban$Q25T[tban$P_ASSIGN3==3& tban$Q11T==0]-1.96*tban$se[tban$P_ASSIGN3==3& tban$Q11T==0],
            tban$Q25T[tban$P_ASSIGN3==2& tban$Q11T==0]-1.96*tban$se[tban$P_ASSIGN3==2& tban$Q11T==0],
            tban$Q25T[tban$P_ASSIGN3==4& tban$Q11T==0]-1.96*tban$se[tban$P_ASSIGN3==4& tban$Q11T==0],0,
            tban$Q25T[tban$P_ASSIGN3==1& tban$Q11T==1]-1.96*tban$se[tban$P_ASSIGN3==1& tban$Q11T==1], 
            tban$Q25T[tban$P_ASSIGN3==3& tban$Q11T==1]-1.96*tban$se[tban$P_ASSIGN3==3& tban$Q11T==1],
            tban$Q25T[tban$P_ASSIGN3==2& tban$Q11T==1]-1.96*tban$se[tban$P_ASSIGN3==2& tban$Q11T==1],
            tban$Q25T[tban$P_ASSIGN3==4& tban$Q11T==1]-1.96*tban$se[tban$P_ASSIGN3==4& tban$Q11T==1])
ci.high <- c(iran$Q21T[iran$P_ASSIGN2==1& iran$Q11T==0]+1.96*iran$se[iran$P_ASSIGN2==1& iran$Q11T==0], 
             iran$Q21T[iran$P_ASSIGN2==4& iran$Q11T==0]+1.96*iran$se[iran$P_ASSIGN2==4& iran$Q11T==0],
             iran$Q21T[iran$P_ASSIGN2==5& iran$Q11T==0]+1.96*iran$se[iran$P_ASSIGN2==5& iran$Q11T==0],0,
             iran$Q21T[iran$P_ASSIGN2==1& iran$Q11T==1]+1.96*iran$se[iran$P_ASSIGN2==1& iran$Q11T==1],
             iran$Q21T[iran$P_ASSIGN2==4& iran$Q11T==1]+1.96*iran$se[iran$P_ASSIGN2==4& iran$Q11T==1],
             iran$Q21T[iran$P_ASSIGN2==5& iran$Q11T==1]+1.96*iran$se[iran$P_ASSIGN2==5& iran$Q11T==1],
             tban$Q25T[tban$P_ASSIGN3==1& tban$Q11T==0]+1.96*tban$se[tban$P_ASSIGN3==1& tban$Q11T==0], 
             tban$Q25T[tban$P_ASSIGN3==3& tban$Q11T==0]+1.96*tban$se[tban$P_ASSIGN3==3& tban$Q11T==0],
             tban$Q25T[tban$P_ASSIGN3==2& tban$Q11T==0]+1.96*tban$se[tban$P_ASSIGN3==2& tban$Q11T==0],
             tban$Q25T[tban$P_ASSIGN3==4& tban$Q11T==0]+1.96*tban$se[tban$P_ASSIGN3==4& tban$Q11T==0],0,
             tban$Q25T[tban$P_ASSIGN3==1& tban$Q11T==1]+1.96*tban$se[tban$P_ASSIGN3==1& tban$Q11T==1], 
             tban$Q25T[tban$P_ASSIGN3==3& tban$Q11T==1]+1.96*tban$se[tban$P_ASSIGN3==3& tban$Q11T==1],
             tban$Q25T[tban$P_ASSIGN3==2& tban$Q11T==1]+1.96*tban$se[tban$P_ASSIGN3==2& tban$Q11T==1],
             tban$Q25T[tban$P_ASSIGN3==4& tban$Q11T==1]+1.96*tban$se[tban$P_ASSIGN3==4& tban$Q11T==1])
plot.labs <- c("Control","Support","Oppose","",
               " Control "," Support "," Oppose ",
               "Control ","Support ","Oppose ","Divided "," ",
               "  Control ","  Support ","  Oppose ","  Divided ")
df.plot <- data.frame(Treatment=plot.labs, estimate=pt.est, 
                      lower=ci.low, upper=ci.high)
df.plot$Treatment <- factor(df.plot$Treatment, levels = plot.labs)

g.i <- ggplot(data=df.plot[c(1:7),], 
              aes(x=Treatment, y=estimate),
              size = 0.5)+
  geom_col(position=position_dodge(1), width=0.75,fill="gray50")+
  geom_errorbar(aes(ymin=lower[1:7],ymax=upper[1:7]),width=0.2)+
  ylim(0,.5)+
  ylab("Respondent Agreement (%)")+
  xlab("Not Confident                                                Confident")+
  ggtitle("Iran Strikes")+
  geom_text(aes(label = round(estimate,3)), vjust = -0.5)+
  theme_bw()

g.tb <- ggplot(data=df.plot[c(8:16),], 
               aes(x=Treatment, y=estimate),
               size = 0.5)+
  geom_col(position=position_dodge(1), width=0.75,fill="gray50")+
  geom_errorbar(aes(ymin=lower[1:9],ymax=upper[1:9]),width=0.2)+ # data is limited to rows 8 to 16, which requires this numbering scheme, instead of 8:16
  ylim(0,.5)+
  ylab("Respondent Agreement (%)")+
  xlab("Not Confident                                                Confident")+
  ggtitle("Transgender Ban")+
  geom_text(aes(label = round(estimate,3)), vjust = -0.5)+
  theme_bw()

grid.arrange(g.i,g.tb,ncol=2)



############################################################################



### Figure 10.4
summary(factor(df$Q25)) # Transgender ban
summary(factor(df$Q11T)) 

df2 <- df
df2 <- df2[df2$Q11 < 77,]

# Transgender Ban treatments
# P_ASSIGN3 == 1    # Control
# P_ASSIGN3 == 3    # Support
# P_ASSIGN3 == 2    # Oppose
# P_ASSIGN3 == 4    # Divided

# Q25T
df2$Q25T <- NA
df2$Q25T[df2$Q25 < 3] <- 1
df2$Q25T[df2$Q25 > 2 & df2$Q25 < 77] <- 0

# Survey designs
w2_design <-
  svydesign(
    id = ~ 1,
    weights = ~ weight2,
    data = df2
  )


# Transgender ban
tban <- svyby(~Q25T,~P_ASSIGN3+party+Q11T, w2_design, svymean, na.rm = TRUE)

# Build out plots
# TBan
pt.est <- c(tban$Q25T[tban$P_ASSIGN3==1& tban$Q11T==0&tban$party==0], 
            tban$Q25T[tban$P_ASSIGN3==3& tban$Q11T==0&tban$party==0],
            tban$Q25T[tban$P_ASSIGN3==2& tban$Q11T==0&tban$party==0],
            tban$Q25T[tban$P_ASSIGN3==4& tban$Q11T==0&tban$party==0],
            0,
            tban$Q25T[tban$P_ASSIGN3==1& tban$Q11T==0&tban$party==2], 
            tban$Q25T[tban$P_ASSIGN3==3& tban$Q11T==0&tban$party==2],
            tban$Q25T[tban$P_ASSIGN3==2& tban$Q11T==0&tban$party==2],
            tban$Q25T[tban$P_ASSIGN3==4& tban$Q11T==0&tban$party==2],
            0,
            tban$Q25T[tban$P_ASSIGN3==1& tban$Q11T==1&tban$party==0], 
            tban$Q25T[tban$P_ASSIGN3==3& tban$Q11T==1&tban$party==0],
            tban$Q25T[tban$P_ASSIGN3==2& tban$Q11T==1&tban$party==0],
            tban$Q25T[tban$P_ASSIGN3==4& tban$Q11T==1&tban$party==0],
            0,
            tban$Q25T[tban$P_ASSIGN3==1& tban$Q11T==1&tban$party==2], 
            tban$Q25T[tban$P_ASSIGN3==3& tban$Q11T==1&tban$party==2],
            tban$Q25T[tban$P_ASSIGN3==2& tban$Q11T==1&tban$party==2],
            tban$Q25T[tban$P_ASSIGN3==4& tban$Q11T==1&tban$party==2])
ci.low <- c(tban$Q25T[tban$P_ASSIGN3==1& tban$Q11T==0&tban$party==0]-1.96*tban$se[tban$P_ASSIGN3==1& tban$Q11T==0&tban$party==0], 
            tban$Q25T[tban$P_ASSIGN3==3& tban$Q11T==0&tban$party==0]-1.96*tban$se[tban$P_ASSIGN3==3& tban$Q11T==0&tban$party==0],
            tban$Q25T[tban$P_ASSIGN3==2& tban$Q11T==0&tban$party==0]-1.96*tban$se[tban$P_ASSIGN3==2& tban$Q11T==0&tban$party==0],
            tban$Q25T[tban$P_ASSIGN3==4& tban$Q11T==0&tban$party==0]-1.96*tban$se[tban$P_ASSIGN3==4& tban$Q11T==0&tban$party==0],0,
            tban$Q25T[tban$P_ASSIGN3==1& tban$Q11T==0&tban$party==2]-1.96*tban$se[tban$P_ASSIGN3==1& tban$Q11T==0&tban$party==2], 
            tban$Q25T[tban$P_ASSIGN3==3& tban$Q11T==0&tban$party==2]-1.96*tban$se[tban$P_ASSIGN3==3& tban$Q11T==0&tban$party==2],
            tban$Q25T[tban$P_ASSIGN3==2& tban$Q11T==0&tban$party==2]-1.96*tban$se[tban$P_ASSIGN3==2& tban$Q11T==0&tban$party==2],
            tban$Q25T[tban$P_ASSIGN3==4& tban$Q11T==0&tban$party==2]-1.96*tban$se[tban$P_ASSIGN3==4& tban$Q11T==0&tban$party==2],0,
            tban$Q25T[tban$P_ASSIGN3==1& tban$Q11T==1&tban$party==0]-1.96*tban$se[tban$P_ASSIGN3==1& tban$Q11T==1&tban$party==0], 
            tban$Q25T[tban$P_ASSIGN3==3& tban$Q11T==1&tban$party==0]-1.96*tban$se[tban$P_ASSIGN3==3& tban$Q11T==1&tban$party==0],
            tban$Q25T[tban$P_ASSIGN3==2& tban$Q11T==1&tban$party==0]-1.96*tban$se[tban$P_ASSIGN3==2& tban$Q11T==1&tban$party==0],
            tban$Q25T[tban$P_ASSIGN3==4& tban$Q11T==1&tban$party==0]-1.96*tban$se[tban$P_ASSIGN3==4& tban$Q11T==1&tban$party==0],0,
            tban$Q25T[tban$P_ASSIGN3==1& tban$Q11T==1&tban$party==2]-1.96*tban$se[tban$P_ASSIGN3==1& tban$Q11T==1&tban$party==2], 
            tban$Q25T[tban$P_ASSIGN3==3& tban$Q11T==1&tban$party==2]-1.96*tban$se[tban$P_ASSIGN3==3& tban$Q11T==1&tban$party==2],
            tban$Q25T[tban$P_ASSIGN3==2& tban$Q11T==1&tban$party==2]-1.96*tban$se[tban$P_ASSIGN3==2& tban$Q11T==1&tban$party==2],
            tban$Q25T[tban$P_ASSIGN3==4& tban$Q11T==1&tban$party==2]-1.96*tban$se[tban$P_ASSIGN3==4& tban$Q11T==1&tban$party==2])
ci.high <- c(tban$Q25T[tban$P_ASSIGN3==1& tban$Q11T==0&tban$party==0]+1.96*tban$se[tban$P_ASSIGN3==1& tban$Q11T==0&tban$party==0], 
             tban$Q25T[tban$P_ASSIGN3==3& tban$Q11T==0&tban$party==0]+1.96*tban$se[tban$P_ASSIGN3==3& tban$Q11T==0&tban$party==0],
             tban$Q25T[tban$P_ASSIGN3==2& tban$Q11T==0&tban$party==0]+1.96*tban$se[tban$P_ASSIGN3==2& tban$Q11T==0&tban$party==0],
             tban$Q25T[tban$P_ASSIGN3==4& tban$Q11T==0&tban$party==0]+1.96*tban$se[tban$P_ASSIGN3==4& tban$Q11T==0&tban$party==0],0,
             tban$Q25T[tban$P_ASSIGN3==1& tban$Q11T==0&tban$party==2]+1.96*tban$se[tban$P_ASSIGN3==1& tban$Q11T==0&tban$party==2], 
             tban$Q25T[tban$P_ASSIGN3==3& tban$Q11T==0&tban$party==2]+1.96*tban$se[tban$P_ASSIGN3==3& tban$Q11T==0&tban$party==2],
             tban$Q25T[tban$P_ASSIGN3==2& tban$Q11T==0&tban$party==2]+1.96*tban$se[tban$P_ASSIGN3==2& tban$Q11T==0&tban$party==2],
             tban$Q25T[tban$P_ASSIGN3==4& tban$Q11T==0&tban$party==2]+1.96*tban$se[tban$P_ASSIGN3==4& tban$Q11T==0&tban$party==2],0,
             tban$Q25T[tban$P_ASSIGN3==1& tban$Q11T==1&tban$party==0]+1.96*tban$se[tban$P_ASSIGN3==1& tban$Q11T==1&tban$party==0], 
             tban$Q25T[tban$P_ASSIGN3==3& tban$Q11T==1&tban$party==0]+1.96*tban$se[tban$P_ASSIGN3==3& tban$Q11T==1&tban$party==0],
             tban$Q25T[tban$P_ASSIGN3==2& tban$Q11T==1&tban$party==0]+1.96*tban$se[tban$P_ASSIGN3==2& tban$Q11T==1&tban$party==0],
             tban$Q25T[tban$P_ASSIGN3==4& tban$Q11T==1&tban$party==0]+1.96*tban$se[tban$P_ASSIGN3==4& tban$Q11T==1&tban$party==0],0,
             tban$Q25T[tban$P_ASSIGN3==1& tban$Q11T==1&tban$party==2]+1.96*tban$se[tban$P_ASSIGN3==1& tban$Q11T==1&tban$party==2], 
             tban$Q25T[tban$P_ASSIGN3==3& tban$Q11T==1&tban$party==2]+1.96*tban$se[tban$P_ASSIGN3==3& tban$Q11T==1&tban$party==2],
             tban$Q25T[tban$P_ASSIGN3==2& tban$Q11T==1&tban$party==2]+1.96*tban$se[tban$P_ASSIGN3==2& tban$Q11T==1&tban$party==2],
             tban$Q25T[tban$P_ASSIGN3==4& tban$Q11T==1&tban$party==2]+1.96*tban$se[tban$P_ASSIGN3==4& tban$Q11T==1&tban$party==2])
plot.labs <- c("Control","Support","Oppose","Divided","",
               "Control ","Support ","Oppose ","Divided "," ",
               " Control"," Support"," Oppose"," Divided","  ",
               " Control "," Support "," Oppose "," Divided ")
party <- c("D","D","D","D","NA","R","R","R","R","NA",
           "D","D","D","D","NA","R","R","R","R")

df.plot <- data.frame(Treatment=plot.labs, estimate=pt.est, 
                      lower=ci.low, upper=ci.high,party=party)
df.plot$Treatment <- factor(df.plot$Treatment, levels = plot.labs)

g.tbp <- ggplot(data=df.plot, 
                aes(x=Treatment, y=estimate),
                size = 0.5)+
  geom_col(position=position_dodge(1), width=0.75,
           fill=c("blue4","blue4",'blue4','blue4','black',"red4","red4","red4","red4","black",
                  "blue4","blue4",'blue4','blue4','black',"red4","red4","red4","red4"))+
  geom_errorbar(aes(ymin=lower,ymax=upper),width=0.2)+
  ylim(0,.8)+
  ylab("Respondent Agreement (%)")+
  xlab("Not Confident                                                                                                             Confident    ")+
  ggtitle("")+
  geom_text(aes(label = round(estimate,3)), vjust = -0.5)

g.tbp



############################################################################



### Figure 10.5
# Treatment DOV_LEAD - 1 Civ 2 Dem 3 Rep

# 1&2 = Agree; 4&5 = Disagree
# Q24A - Political leaders had a good plan
# Q24B - Military leaders implemented as well as they could
# Q24C - Military leaders gave good advice about strategy
# Q24D - Pol leaders listened to mil leaders as much as they should
# Q24E - Pol leaders did a good job integrating DIME in IRQ and AFG
# Q24F - Mil leaders did a good job integrating ops with Dip Econ in I&A

df2 <- df
df2 <- df2[df2$Q11 < 77,]

summary(factor(df2$DOV_LEAD))
# Panel 1
summary(factor(df2$Q24A))
summary(factor(df2$Q24B))

dfp1 <- df2[df2$Q24A < 6 & df2$Q24B < 6,]

dfp1$Q24AT <- 0
dfp1$Q24AT[dfp1$Q24A < 3] <- 1

dfp1$Q24BT <- 0
dfp1$Q24BT[dfp1$Q24B < 3] <- 1

# Panel 2
summary(factor(df2$Q24D))
summary(factor(df2$Q24C))

dfp2 <- df2[df2$Q24D < 6 & df2$Q24C < 6,]

dfp2$Q24DT <- 0
dfp2$Q24DT[dfp2$Q24D < 3] <- 1

dfp2$Q24CT <- 0
dfp2$Q24CT[dfp2$Q24C < 3] <- 1

# Panel 3
summary(factor(df2$Q24E))
summary(factor(df2$Q24F))

dfp3 <- df2[df2$Q24E < 6 & df2$Q24F < 6,]

dfp3$Q24ET <- 0
dfp3$Q24ET[dfp3$Q24E < 3] <- 1

dfp3$Q24FT <- 0
dfp3$Q24FT[dfp3$Q24F < 3] <- 1


# Survey designs
w2_design_p1 <-
  svydesign(
    id = ~ 1,
    weights = ~ weight2,
    data = dfp1
  )

w2_design_p2 <-
  svydesign(
    id = ~ 1,
    weights = ~ weight2,
    data = dfp2
  )

w2_design_p3 <-
  svydesign(
    id = ~ 1,
    weights = ~ weight2,
    data = dfp3
  )

# Transgender ban
cp1 <- svyby(~Q24AT,~DOV_LEAD+party, w2_design_p1, svymean, na.rm = TRUE)
mp1 <- svyby(~Q24BT,~party, w2_design_p1, svymean, na.rm = TRUE)
mp1$DOV_LEAD <- 4
colnames(cp1)[3] <- "agreement"
colnames(mp1)[2] <- "agreement"
p1 <- rbind(cp1,mp1)

cp2 <- svyby(~Q24DT,~DOV_LEAD+party, w2_design_p2, svymean, na.rm = TRUE)
mp2 <- svyby(~Q24CT,~party, w2_design_p2, svymean, na.rm = TRUE)
mp2$DOV_LEAD <- 4
colnames(cp2)[3] <- "agreement"
colnames(mp2)[2] <- "agreement"
p2 <- rbind(cp2,mp2)

cp3 <- svyby(~Q24ET,~DOV_LEAD+party, w2_design_p3, svymean, na.rm = TRUE)
mp3 <- svyby(~Q24FT,~party, w2_design_p3, svymean, na.rm = TRUE)
mp3$DOV_LEAD <- 4
colnames(cp3)[3] <- "agreement"
colnames(mp3)[2] <- "agreement"
p3 <- rbind(cp3,mp3)


# Build out plots
# plot 1
pt.est1 <- c(p1$agreement[p1$party==0&p1$DOV_LEAD==1], 
             p1$agreement[p1$party==0&p1$DOV_LEAD==2],
             p1$agreement[p1$party==0&p1$DOV_LEAD==3],
             p1$agreement[p1$party==0&p1$DOV_LEAD==4],
             0,
             p1$agreement[p1$party==1&p1$DOV_LEAD==1], 
             p1$agreement[p1$party==1&p1$DOV_LEAD==2],
             p1$agreement[p1$party==1&p1$DOV_LEAD==3],
             p1$agreement[p1$party==1&p1$DOV_LEAD==4],
             0,
             p1$agreement[p1$party==2&p1$DOV_LEAD==1], 
             p1$agreement[p1$party==2&p1$DOV_LEAD==2],
             p1$agreement[p1$party==2&p1$DOV_LEAD==3],
             p1$agreement[p1$party==2&p1$DOV_LEAD==4])
ci.low1 <- c(p1$agreement[p1$party==0&p1$DOV_LEAD==1]-1.96*p1$se[p1$party==0&p1$DOV_LEAD==1], 
             p1$agreement[p1$party==0&p1$DOV_LEAD==2]-1.96*p1$se[p1$party==0&p1$DOV_LEAD==2],
             p1$agreement[p1$party==0&p1$DOV_LEAD==3]-1.96*p1$se[p1$party==0&p1$DOV_LEAD==3],
             p1$agreement[p1$party==0&p1$DOV_LEAD==4]-1.96*p1$se[p1$party==0&p1$DOV_LEAD==4],
             0,
             p1$agreement[p1$party==1&p1$DOV_LEAD==1]-1.96*p1$se[p1$party==1&p1$DOV_LEAD==1], 
             p1$agreement[p1$party==1&p1$DOV_LEAD==2]-1.96*p1$se[p1$party==1&p1$DOV_LEAD==2],
             p1$agreement[p1$party==1&p1$DOV_LEAD==3]-1.96*p1$se[p1$party==1&p1$DOV_LEAD==3],
             p1$agreement[p1$party==1&p1$DOV_LEAD==4]-1.96*p1$se[p1$party==1&p1$DOV_LEAD==4],
             0,
             p1$agreement[p1$party==2&p1$DOV_LEAD==1]-1.96*p1$se[p1$party==2&p1$DOV_LEAD==1], 
             p1$agreement[p1$party==2&p1$DOV_LEAD==2]-1.96*p1$se[p1$party==2&p1$DOV_LEAD==2],
             p1$agreement[p1$party==2&p1$DOV_LEAD==3]-1.96*p1$se[p1$party==2&p1$DOV_LEAD==3],
             p1$agreement[p1$party==2&p1$DOV_LEAD==4]-1.96*p1$se[p1$party==2&p1$DOV_LEAD==4])
ci.high1 <- c(p1$agreement[p1$party==0&p1$DOV_LEAD==1]+1.96*p1$se[p1$party==0&p1$DOV_LEAD==1], 
              p1$agreement[p1$party==0&p1$DOV_LEAD==2]+1.96*p1$se[p1$party==0&p1$DOV_LEAD==2],
              p1$agreement[p1$party==0&p1$DOV_LEAD==3]+1.96*p1$se[p1$party==0&p1$DOV_LEAD==3],
              p1$agreement[p1$party==0&p1$DOV_LEAD==4]+1.96*p1$se[p1$party==0&p1$DOV_LEAD==4],
              0,
              p1$agreement[p1$party==1&p1$DOV_LEAD==1]+1.96*p1$se[p1$party==1&p1$DOV_LEAD==1], 
              p1$agreement[p1$party==1&p1$DOV_LEAD==2]+1.96*p1$se[p1$party==1&p1$DOV_LEAD==2],
              p1$agreement[p1$party==1&p1$DOV_LEAD==3]+1.96*p1$se[p1$party==1&p1$DOV_LEAD==3],
              p1$agreement[p1$party==1&p1$DOV_LEAD==4]+1.96*p1$se[p1$party==1&p1$DOV_LEAD==4],
              0,
              p1$agreement[p1$party==2&p1$DOV_LEAD==1]+1.96*p1$se[p1$party==2&p1$DOV_LEAD==1], 
              p1$agreement[p1$party==2&p1$DOV_LEAD==2]+1.96*p1$se[p1$party==2&p1$DOV_LEAD==2],
              p1$agreement[p1$party==2&p1$DOV_LEAD==3]+1.96*p1$se[p1$party==2&p1$DOV_LEAD==3],
              p1$agreement[p1$party==2&p1$DOV_LEAD==4]+1.96*p1$se[p1$party==2&p1$DOV_LEAD==4])
plot.labs1 <- c("Civilian","Democratic","Republican","Military","",
                "Civilian ","Democratic ","Republican ","Military "," ",
                " Civilian"," Democratic"," Republican"," Military")

df.plot1 <- data.frame(Treatment=plot.labs1, estimate=pt.est1, 
                       lower=ci.low1, upper=ci.high1)
df.plot1$Treatment <- factor(df.plot1$Treatment, levels = plot.labs1)

g.p1 <- ggplot(data=df.plot1, 
               aes(x=Treatment, y=estimate),
               size = 0.5)+
  geom_col(position=position_dodge(1), width=0.75)+
  geom_errorbar(aes(ymin=lower,ymax=upper),width=0.2)+
  ylim(0,.65)+
  ylab("Respondent Agreement (%)")+
  xlab("\nDemocratic Respondents                                                                            Independent Respondents                                                                           Republican Respondents")+
  ggtitle("Political leaders had a good plan/Military leaders implemented plan well")+
  geom_text(aes(label = round(estimate,3)), vjust = -0.5)

# plot 2
pt.est2 <- c(p2$agreement[p2$party==0&p2$DOV_LEAD==1], 
             p2$agreement[p2$party==0&p2$DOV_LEAD==2],
             p2$agreement[p2$party==0&p2$DOV_LEAD==3],
             p2$agreement[p2$party==0&p2$DOV_LEAD==4],
             0,
             p2$agreement[p2$party==1&p2$DOV_LEAD==1], 
             p2$agreement[p2$party==1&p2$DOV_LEAD==2],
             p2$agreement[p2$party==1&p2$DOV_LEAD==3],
             p2$agreement[p2$party==1&p2$DOV_LEAD==4],
             0,
             p2$agreement[p2$party==2&p2$DOV_LEAD==1], 
             p2$agreement[p2$party==2&p2$DOV_LEAD==2],
             p2$agreement[p2$party==2&p2$DOV_LEAD==3],
             p2$agreement[p2$party==2&p2$DOV_LEAD==4])
ci.low2 <- c(p2$agreement[p2$party==0&p2$DOV_LEAD==1]-1.96*p2$se[p2$party==0&p2$DOV_LEAD==1], 
             p2$agreement[p2$party==0&p2$DOV_LEAD==2]-1.96*p2$se[p2$party==0&p2$DOV_LEAD==2],
             p2$agreement[p2$party==0&p2$DOV_LEAD==3]-1.96*p2$se[p2$party==0&p2$DOV_LEAD==3],
             p2$agreement[p2$party==0&p2$DOV_LEAD==4]-1.96*p2$se[p2$party==0&p2$DOV_LEAD==4],
             0,
             p2$agreement[p2$party==1&p2$DOV_LEAD==1]-1.96*p2$se[p2$party==1&p2$DOV_LEAD==1], 
             p2$agreement[p2$party==1&p2$DOV_LEAD==2]-1.96*p2$se[p2$party==1&p2$DOV_LEAD==2],
             p2$agreement[p2$party==1&p2$DOV_LEAD==3]-1.96*p2$se[p2$party==1&p2$DOV_LEAD==3],
             p2$agreement[p2$party==1&p2$DOV_LEAD==4]-1.96*p2$se[p2$party==1&p2$DOV_LEAD==4],
             0,
             p2$agreement[p2$party==2&p2$DOV_LEAD==1]-1.96*p2$se[p2$party==2&p2$DOV_LEAD==1], 
             p2$agreement[p2$party==2&p2$DOV_LEAD==2]-1.96*p2$se[p2$party==2&p2$DOV_LEAD==2],
             p2$agreement[p2$party==2&p2$DOV_LEAD==3]-1.96*p2$se[p2$party==2&p2$DOV_LEAD==3],
             p2$agreement[p2$party==2&p2$DOV_LEAD==4]-1.96*p2$se[p2$party==2&p2$DOV_LEAD==4])
ci.high2 <- c(p2$agreement[p2$party==0&p2$DOV_LEAD==1]+1.96*p2$se[p2$party==0&p2$DOV_LEAD==1], 
              p2$agreement[p2$party==0&p2$DOV_LEAD==2]+1.96*p2$se[p2$party==0&p2$DOV_LEAD==2],
              p2$agreement[p2$party==0&p2$DOV_LEAD==3]+1.96*p2$se[p2$party==0&p2$DOV_LEAD==3],
              p2$agreement[p2$party==0&p2$DOV_LEAD==4]+1.96*p2$se[p2$party==0&p2$DOV_LEAD==4],
              0,
              p2$agreement[p2$party==1&p2$DOV_LEAD==1]+1.96*p2$se[p2$party==1&p2$DOV_LEAD==1], 
              p2$agreement[p2$party==1&p2$DOV_LEAD==2]+1.96*p2$se[p2$party==1&p2$DOV_LEAD==2],
              p2$agreement[p2$party==1&p2$DOV_LEAD==3]+1.96*p2$se[p2$party==1&p2$DOV_LEAD==3],
              p2$agreement[p2$party==1&p2$DOV_LEAD==4]+1.96*p2$se[p2$party==1&p2$DOV_LEAD==4],
              0,
              p2$agreement[p2$party==2&p2$DOV_LEAD==1]+1.96*p2$se[p2$party==2&p2$DOV_LEAD==1], 
              p2$agreement[p2$party==2&p2$DOV_LEAD==2]+1.96*p2$se[p2$party==2&p2$DOV_LEAD==2],
              p2$agreement[p2$party==2&p2$DOV_LEAD==3]+1.96*p2$se[p2$party==2&p2$DOV_LEAD==3],
              p2$agreement[p2$party==2&p2$DOV_LEAD==4]+1.96*p2$se[p2$party==2&p2$DOV_LEAD==4])
plot.labs2 <- c("Civilian","Democratic","Republican","Military","",
                "Civilian ","Democratic ","Republican ","Military "," ",
                " Civilian"," Democratic"," Republican"," Military")

df.plot2 <- data.frame(Treatment=plot.labs2, estimate=pt.est2, 
                       lower=ci.low2, upper=ci.high2)
df.plot2$Treatment <- factor(df.plot2$Treatment, levels = plot.labs2)

g.p2 <- ggplot(data=df.plot2, 
               aes(x=Treatment, y=estimate),
               size = 0.5)+
  geom_col(position=position_dodge(1), width=0.75)+
  geom_errorbar(aes(ymin=lower,ymax=upper),width=0.2)+
  ylim(0,.65)+
  ylab("Respondent Agreement (%)")+
  xlab("\nDemocratic Respondents                                                                            Independent Respondents                                                                           Republican Respondents")+
  ggtitle("Political leaders listened to military advice/Military leaders gave good advice")+
  geom_text(aes(label = round(estimate,3)), vjust = -0.5)

# plot 3
pt.est3 <- c(p3$agreement[p3$party==0&p3$DOV_LEAD==1], 
             p3$agreement[p3$party==0&p3$DOV_LEAD==2],
             p3$agreement[p3$party==0&p3$DOV_LEAD==3],
             p3$agreement[p3$party==0&p3$DOV_LEAD==4],
             0,
             p3$agreement[p3$party==1&p3$DOV_LEAD==1], 
             p3$agreement[p3$party==1&p3$DOV_LEAD==2],
             p3$agreement[p3$party==1&p3$DOV_LEAD==3],
             p3$agreement[p3$party==1&p3$DOV_LEAD==4],
             0,
             p3$agreement[p3$party==2&p3$DOV_LEAD==1], 
             p3$agreement[p3$party==2&p3$DOV_LEAD==2],
             p3$agreement[p3$party==2&p3$DOV_LEAD==3],
             p3$agreement[p3$party==2&p3$DOV_LEAD==4])
ci.low3 <- c(p3$agreement[p3$party==0&p3$DOV_LEAD==1]-1.96*p3$se[p3$party==0&p3$DOV_LEAD==1], 
             p3$agreement[p3$party==0&p3$DOV_LEAD==2]-1.96*p3$se[p3$party==0&p3$DOV_LEAD==2],
             p3$agreement[p3$party==0&p3$DOV_LEAD==3]-1.96*p3$se[p3$party==0&p3$DOV_LEAD==3],
             p3$agreement[p3$party==0&p3$DOV_LEAD==4]-1.96*p3$se[p3$party==0&p3$DOV_LEAD==4],
             0,
             p3$agreement[p3$party==1&p3$DOV_LEAD==1]-1.96*p3$se[p3$party==1&p3$DOV_LEAD==1], 
             p3$agreement[p3$party==1&p3$DOV_LEAD==2]-1.96*p3$se[p3$party==1&p3$DOV_LEAD==2],
             p3$agreement[p3$party==1&p3$DOV_LEAD==3]-1.96*p3$se[p3$party==1&p3$DOV_LEAD==3],
             p3$agreement[p3$party==1&p3$DOV_LEAD==4]-1.96*p3$se[p3$party==1&p3$DOV_LEAD==4],
             0,
             p3$agreement[p3$party==2&p3$DOV_LEAD==1]-1.96*p3$se[p3$party==2&p3$DOV_LEAD==1], 
             p3$agreement[p3$party==2&p3$DOV_LEAD==2]-1.96*p3$se[p3$party==2&p3$DOV_LEAD==2],
             p3$agreement[p3$party==2&p3$DOV_LEAD==3]-1.96*p3$se[p3$party==2&p3$DOV_LEAD==3],
             p3$agreement[p3$party==2&p3$DOV_LEAD==4]-1.96*p3$se[p3$party==2&p3$DOV_LEAD==4])
ci.high3 <- c(p3$agreement[p3$party==0&p3$DOV_LEAD==1]+1.96*p3$se[p3$party==0&p3$DOV_LEAD==1], 
              p3$agreement[p3$party==0&p3$DOV_LEAD==2]+1.96*p3$se[p3$party==0&p3$DOV_LEAD==2],
              p3$agreement[p3$party==0&p3$DOV_LEAD==3]+1.96*p3$se[p3$party==0&p3$DOV_LEAD==3],
              p3$agreement[p3$party==0&p3$DOV_LEAD==4]+1.96*p3$se[p3$party==0&p3$DOV_LEAD==4],
              0,
              p3$agreement[p3$party==1&p3$DOV_LEAD==1]+1.96*p3$se[p3$party==1&p3$DOV_LEAD==1], 
              p3$agreement[p3$party==1&p3$DOV_LEAD==2]+1.96*p3$se[p3$party==1&p3$DOV_LEAD==2],
              p3$agreement[p3$party==1&p3$DOV_LEAD==3]+1.96*p3$se[p3$party==1&p3$DOV_LEAD==3],
              p3$agreement[p3$party==1&p3$DOV_LEAD==4]+1.96*p3$se[p3$party==1&p3$DOV_LEAD==4],
              0,
              p3$agreement[p3$party==2&p3$DOV_LEAD==1]+1.96*p3$se[p3$party==2&p3$DOV_LEAD==1], 
              p3$agreement[p3$party==2&p3$DOV_LEAD==2]+1.96*p3$se[p3$party==2&p3$DOV_LEAD==2],
              p3$agreement[p3$party==2&p3$DOV_LEAD==3]+1.96*p3$se[p3$party==2&p3$DOV_LEAD==3],
              p3$agreement[p3$party==2&p3$DOV_LEAD==4]+1.96*p3$se[p3$party==2&p3$DOV_LEAD==4])
plot.labs3 <- c("Civilian","Democratic","Republican","Military","",
                "Civilian ","Democratic ","Republican ","Military "," ",
                " Civilian"," Democratic"," Republican"," Military")

df.plot3 <- data.frame(Treatment=plot.labs3, estimate=pt.est3, 
                       lower=ci.low3, upper=ci.high3)
df.plot3$Treatment <- factor(df.plot3$Treatment, levels = plot.labs3)

g.p3 <- ggplot(data=df.plot3, 
               aes(x=Treatment, y=estimate),
               size = 0.5)+
  geom_col(position=position_dodge(1), width=0.75)+
  geom_errorbar(aes(ymin=lower,ymax=upper),width=0.2)+
  ylim(0,.65)+
  ylab("Respondent Agreement (%)")+
  xlab("\nDemocratic Respondents                                                                            Independent Respondents                                                                           Republican Respondents")+
  ggtitle("Political/Military leaders did a good job integrating military, diplomatic, and economic efforts")+
  geom_text(aes(label = round(estimate,3)), vjust = -0.5)

grid.arrange(g.p1,g.p2,g.p3,ncol=1,nrow=3)

